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ABSTRACT 

We use 3D hydro dynamical models to investigate the effects of massive star feedback 
from winds and supernovae on inhomogeneous molecular material left over from the 
formation of a massive stellar cluster. We simulate the interaction of the mechanical 
energy input from a cluster with 3 0-stars into a giant molecular cloud (GMC) clump 
containing 3240 of molecular material within a 4 pc radius. The cluster wind blows 
out of the molecular clump along low-density channels, into which denser clump mate- 
rial is entrained. We find that the densest molecular regions are surprisingly resistant 
to ablation by the cluster wind, in part due to shielding by other dense regions closer 
to the cluster. Nonetheless, molecular material is gradually removed by the cluster 
wind during which mass-loading factors in excess of several 100 are obtained. Because 
the clump is very porous, 60 — 75 per cent of the injected wind energy escapes the 
simulation domain, with the difference being radiated. After 4.4 Myr, the massive stars 
in our simulation begin to explode as supernovae. The highly structured environment 
into which the SN energy is released allows even weaker coupling to the remaining 
dense material and practically all of the SN energy reaches the wider environment. 
The molecular material is almost completely dispersed and destroyed after 6 Myr. 
The escape fraction of ionizing radiation is estimated to be about 50 per cent during 
the first 4 Myr of the cluster's life. A similar model with a larger and more massive 
GMC clump reveals the same general picture, though more time is needed for it to be 
destroyed. 
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1 INTRODUCTION 

Massive stars have a profound effect on their natal environ- 
ment creating wind-blown shells, cavities and HIT regions. 
Their winds and supernovae (SNe) chemically enrich the in- 
terstellar medium (ISM) and also help to sustain turbulence 
within it. Massive stars embedded within molecular clouds 
likely inhibit further star formation as their winds and ioniz- 
ing radiation disperse and destroy the remaining molecular 
gas, though in some circumstances massive stars may also 



trigger new star formation (Koenig et al. 2012|) and new 
cluster formation (Beuther et al. 2008; Gray & Scannapieco 



2011). The removal of molecular material is also crucial to 



the question of cluster dissolution (Portegies Zwart et al. 


2010 Pfalzner||2011| 


Pelupessy & Portegies Zwart|2012 


). 



Stellar feedback is also recognized as having significant 
influence on galactic and extragalactic scales. For instance, 
without strong stellar feedback, cosmological models predict 
around 10 times the stellar mass found in real galaxies (e.g. 
Cole et al.||2000||Keres et al.pOOQ ). Feedback from massive 
stars can also drive galactic winds from starburst galaxies 



(e.g. [Axon Taylor||1978| [Bland Tulley|[T988l [Heckman 



et al.[[2000 Adelberger et al. 2003), and appears to be re- 
sponsible for the low star formation efficiency of galaxies 
with dark matter haloes somewhat less than the halo mass 
of the Milky Way ( |Guo et al.|2010 ). Ionizing radiation from 
massive stars is also important for the ionization of galaxies 
( Reynolds [1984 ) and the reionization of the early universe 



(Fan et al.[[2006). 



In order to understand the effects of massive stars on 
galactic scales, however, we must first understand their im- 
pact on their local (cluster) environment. The extent to 
which a cloud is affected by stellar feedback is clearly de- 
pendent on a number of parameters, including the mass of 
the cloud and the stellar cluster, the structure of the cloud, 
the position of the cluster relative to the cloud, and the age 
of the system. However, the degree to which stellar feedback 
processes (stellar winds, SNe, and ionizing radiation) cou- 
ple to the clumpy, inhomogeneous molecular clouds which 
initially surround a massive stellar cluster is exceedingly ill- 
determined, and the dominant feedback process is still to 



be settled (e.g. Yorke et al. 1989 Draine k Woods 1991 
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Matzner" 2002 ; Wang et aL][20Tol [Lopez et ar]|20T7] |Pelle-| 
grini et al. 2011J . 



That stellar winds play a significant role in stellar clus- 
ter feedback is apparent from the fact that many young (pre- 
supernova) massive star forming regions contain diffuse X- 
ray emission - the hot gas responsible can only have been 
created by winds. Observations reveal that the surround- 
ing cold molecular material can sometimes confine this hot 
gas (e.g. [Townsley et al.^'2006), but around other clusters 
the cold clouds appear to be shaped and removed by the 
hot gas. In fact, there are now several lines of evidence that 
indicate that the hot X-ray emitting gas often escapes or 
leaks out of the local cluster environment, instead of being 
bottled up inside a swept-up shell. Firstly, the diffuse X-ray 
emission in Ml 7 and the Rosette nebula reveals that only 
a small proportion of the cluster wind energy is radiated in 



the X-ray, and Townsley et al. ( 2003 ) conclude that most of 



the hot gas must flow without cooling into the wider ISM. 
This picture is supported by analysis of the Omega nebula, 
the Arches cluster, and NGC3603, all of which contain an 
amount of hot gas equal to only 1 per cent of the wind ma- 
terial from the 0-stars over the age of the cluster. Secondly, 
direct evidence for outflowing gas comes from observations 
of M17 and particularly RCW49 which show stellar bow 
shocks around 0-stars outside of the central clusters. These 
indicate large scale gas outflows away from the stell ar clus- 
ter with velocities of at least a few hundred kms~^ (Povich 
et al.||2008t 



Stellar wind feedback into an inhomogeneous environ- 
ment has been considered by Tan & McKee ( 2001| ) and 
more recently by|Harper-Clark & Murray (2009). The latter 



postulate that the non-uniform surrounding medium causes 
gaps in the swept-up shell surrounding the wind-blown bub- 
ble where some of the high-pressure gas in the bubble inte- 
rior can leak out. This scenario has received backing from 



Lopez et al. (2011) who conclude that such leakage may be 
occuring within 30 Doradus. Lope z et al.| |2011") also con- 
clude that direct stellar radiation pressure dominates the 
interior dynamics, but this claim has proved far more con- 
troversial, and other works argue in favour of the thermal 
pressure of hot X-ray emitting plasma shaping the large- 
scale structure and dynamics in 30 Doradus (Pellegrini et al. 
2011 Povich 2012). One must recognize that there are sig- 



nificant uncertainties in determining the hot gas pressure in 
regions like 30 Doradus, where it is unclear whether the X- 
ray emitting gas should be treated as a single large bubble 
or as multiple smaller bubbles with distinct identities. 

The possibility that the pressure exerted by stellar radi- 
ation may be dynamically important in massive young stel- 
lar clusters has received much attention in recent years, with 
Krumho lz k Matznerl (|2009l, ^Fall et al. (2010) and Murray 
et al. (2010) all arguing that radiation pressure is the dom- 
inant feedback mechanism. However, these works disagree 
on the net momentum coupling between the radiation field 
and the gas, partially because this depends on the degree of 
inhomogeneity of the gas and the effect that this has on the 
radiation field (e.g. Krumholz Thompson|2012 ). Comple- 
mentary work on the ionized gas pressure has shown that 
ionization feedback into a highly inhomogeneous medium 
is not very effective at high cluster masses (Dale & Bonnell 
2011 ), but becomes more so at lower masses ( |Dale, Ercolano, 
&; Bonnell|2'0T2t . 



Given these competing processes, our aim in this work 
is to examine the extent to which the mechanical energy in- 
put from a cluster of massive stars is confined by and shapes 
the local environment. In this initial investigation we focus 
solely on feedback due to stellar winds and supernovae, and 
defer investigations of other forms of feedback (radiation 
pressure, photoionization etc.) to future works. We conduct 
our investigation through 3D hydrodynamic simulations of 
this interaction. In Section|2]we describe the numerical mod- 
els and initial conditions used in this work. We present and 
discuss results from the simulations in Section |3] Section |4] 
summarizes and concludes this work. 



2 SIMULATIONS OF STELLAR FEEDBACK 

We use an MPI-parallelized numerical scheme to solve the 
Euler equations of hydrodynamics using a lagrangian for- 
mulation and a remap onto the original grid. A piecewise 
parabolic interpolation and characteristic tracing is used to 
obtain the time-averaged fluid variables at each zone inter- 
face. The code then solves a Riemann problem to determine 
the time-averaged fluxes, and then solves the equations of 
hydrodynamics: 



dp 
dt 



+ V • (pu) 



dpu 

dpe 



+ V • {puu + P) 
+ V • [{pe + P)u] 



0, 
0, 

nT — n^A, 



(1) 

(2) 
(3) 



where e — u^/2+e/p is the total specific energy, p is the mass 
density, e is the internal energy density, P is the pressure, 
and T is the temperature. We adopt an ideal gas equation 
of state, e = P/{l — 1), and solar abundances. 

The net heating/cooling rate per unit volume is param- 
eterized as e = nP — n^A, where n — p/mu^ and P and A 
are heating and cooling coefficients which are assumed to 
depend only on temperature. In the ISM, P decreases with 
increasing density as the starlight, soft X-ray, and cosmic 
ray flux are attenuated by the high column densities as- 
sociated with dense clouds. Because the exact form of the 
attenuation depends on details which remain uncertain (for 
instance the size and abundance of PAHs) , the heating rate 
at T < 10^ K is similarly uncertain. In this work we as- 
sume that P = 10"^^ erg s"^ (independent of p or T). The 
low-temperature (T < 10^ K) cooling was then adjusted 
to give 3 thermally stable phases at thermal pressures be- 
tween 2000 — 6000 cm~^ K, as required by observations. 
These stable phases, at temperatures ^ lOK, ~ 150 K and 
^ 8500 K, correspond to the molecular, atomic and warm 
neutral/ionized phases, respectively. The cooling curve and 
phase diagram are shown in iPittard (2011). 



The simulation uses a temperature-dependent average 
particle mass, p. In the molecular phase p = 2.36, reducing 
to 0.61 in ionized gas. The value of p is determined from a 
look-up table of values oip/ p ( |Sutherland|2010 ) . A temper- 
ature independent value of 7, the ratio of specific heats, is 
used, and we set 7 = 5/3. 

Simulations are performed on a 512^ grid with free out- 
fiow boundary conditions. A number of advected scalars are 
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Table 1. Wind properties of the three stars in the cluster as they 
evolve. 
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included to trace the different origins of the gas - stellar 
wind/SN, GMC clump, and surrounding ambient ISM. 



2.1 Stellar Evolution 

We assume that the stellar wind feedback into the surround- 
ing GMC clump is dominated by 3 O stars with initial 
masses of 35 M©, 32 M© and 28 M©, and Main Sequence 
(MS) mass-loss rates of 5, 2.5 and 1.5x 10~^ M© yr~^ re- 
spectively. For each star the MS wind terminal velocity is 
assumed to be 2000 km s~^. Three evolutionary phases are 
considered for each star. The MS phase of the 35 M© star 
lasts for 4 Myr, after which the star becomes a Red Super- 
giant (RSG) and blows a dense, very slow wind character- 
ized by M = lO"'^ Mq yr~^ and ^;oo = 50 km s~^ This phase 
lasts for 0.1 Myr after which the star enters the Wolf-Rayet 
(WR) stage, and a fast, high momentum wind is blown 
(M = 2 x IQ-^Moyr-^ and ^;oo = 2000 km s"^). The WR 
phase lasts for a further 0.3 Myr, after which the star un- 
dergoes a supernova explosion. At this point the star's wind 
is switched off, and 10^^ ergs of thermal energy is imparted 
to the simulation along with 10 M© of ejecta. The other two 
stars remain on the MS throughout this entire evolution- 
ary period. The 32 M© star evolves onto the RSG branch 
0.1 Myr after the 35 M© star explodes. The wind values and 
lifetimes used are summarized in Table [T] and are intended 
to be representative of stars of these masses. 



2.2 Initial Conditions 

The cluster wind blows into a turbulent and inhomoge- 
neous GMC clump, whose structure is based on the work 
of [Vazquez- Semadeni et al. ( 2008) of turbulent and clumpy 
molecular clouds (specifically model Ms24J6). This model 
has a nominal Mach number of 15, is isothermal, and has 
no magnetic field. We scale these results to create a GMC 
clump of radius 4pc and mass 3240 M© (SimA). This gives 



of density 3.33 x 10"^^gcm"^ (nn ^ r 
temperature 8000 K. 

The hydrodynamic grid covers a cubic region of ±16 pc 
extent centered on the GMC clump. The cluster wind is in- 
jected as purely thermal energy within a radius of 0.375 pc (6 
cells). The densest regions cool to 1 K, which is the imposed 
temperature fioor. 

We also explore another model where the GMC clump 
radius is 5pc and where the density and pressure of the 
GMC clump and its surroundings are twice as great (SimB). 
This produces a clump mass of 10,500 Mq (due to the inho- 
mogeneity of the clump it is not quite 3.9 x as massive as 
the clump in SimA). Unless otherwise noted, all results are 
for SimA. 



2.3 Neglected Processes and Simplifications 

This work is the first step in examining the feedback from a 
stellar cluster into surrounding molecular material left over 
from its formation. As such we necessarily make many sim- 
plifications and approximations. 

Our simulations do not include gravity, thermal conduc- 
tion, magnetic fields, radiation pressure, photoionization, 
dust, cosmic rays or radiative losses within the stellar cluster. 
We estimate from the parameters noted in Table ^ that the 
winds have comparable momentum to the radiation fields 
of the stars (both emit momentum of roughly 10^"^ gcms~^ 
over the lifetime of the clusteiQ . The key factor in determin- 
ing the relative importance of the winds and radiation fields 
to the dynamics of the gas is the strength of the coupling 
of the radiation field to the gas. If the degree of wind and 
photon leakage out of the cluster is comparable, then radia- 
tion pressure will provide just an order-unity enhancement 
to the wind pressure ( Krumholz &: Matzner|2009 ). We note 
that our neglect of direct radiation pressure and gravity off- 
set each other to some degree. However, the self- gravity of 
some of the densest structures in our cluster may be impor- 
tant on the timescales that we consider. 

Inclusion of thermal conduction and photoevaporation 
should speed up the destruction of molecular material. Our 
simulations reveal that the densest cloud fragments have 
10~^^ gcm~^ and radii of about 0.1 pc. At a distance of 



p ^ lU ^ g cm 

3 pc from a stellar cluster emitting 10^° ionizing photons per 
second, we estimate that the mass- loss rate from photoevap- 
oration(see, e.g., Pittard|2007 for the relevant equations) is 



4 x 10^ 



an average density of?^8xl0~ g cm~ , or a molecular hy- 
drogen number density nB.2 ~ 250 cm""^. The clump initially 
has a uniform temperature of about 10 K, and is in rough 



'gs giving a lifetime of about 1.5 Myr. This 
is comparable to the mass-loss rate and lifetime due to hy- 
drodynamic ablation. So we should expect slightly quicker 
destruction of molecular material than occurs in our sim- 
ulations through ablation alone. Having said this, photoe- 
vaporation may be suppressed in regions where the ram or 
thermal pressure of the surrounding medium is greater than 
the pressure of the evaporating flow (Dyson 1994). Given 
these considerations, we believe that our simulations should 



^ Note that in an optically thick medium where each photon is 
absorbed and reemitted multiple times, the momentum deposited 
is limited by the energy rather than by the momentum of the 
radiation field. 
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be reasonably representative even with our neglect of photo- 
evaporation. Neglecting thermal conduction may also make 
little difference to our results given that it will be confined 
to regions where there is a temperature gradient parallel to 
the magnetic field. 

We also note that dust can dramatically affect the cool- 
ing within hot bubbles if it can be continuously replenished, 
perhaps by the evaporation/destruction of dense clumps 
( Everett Churchwell||2010|) wh ich will also mass- load the 
bubble (see, e.g., [Pittard et al.| |2Q01a|bt . The presence of 
dust will also affect the photoionization rate throughout the 
cluster. Clearly the effects of dust warrant study in future 
work. 

Radiative losses within the cluster can significantly 
change the energy flux into the surrounding environment 
( Silich et al.||2004| ). However, this is not a significant effect 
for the cluster parameters we have considered, and is likely 
important only for the most massive clusters and super star 
clusters. 

A final consideration concerns the initial conditions. In 
this first work we have used the results of simulations of a 
turbulent ISM, truncated at a radius of 4 or 5 pc around the 
stellar cluster. In future we will attempt to construct mod- 
els which are more self-consistent, by investigating feedback 
into structures formed from colliding flows, clouds and fila- 
ments. 



3 RESULTS 

In this section we present our results. We examine the ini- 
tial blowout (Sec. 3.1), then the evolution up to t = 4Myr 
we focus on the feedback during the 



3.3 



(Sec. [3^. In Sec. 
evolution of the stars through their various evolved stages 
(RSG and WR) and their subsequent explosions as SNe. We 
examine the mass and energy fluxes into the wider surround- 
ings in Sec. |3.4| and the evolution of the column density in 
Sec. |3.5| A comparison between SimA and SimB, plus the 
evolution of the molecular mass is made in Sec. [ 



3.1 Initial Blowout 

The cluster wind creates a high pressure and high temper- 
ature bubble within the GMC clump which expands most 
rapidly into regions of lower density. The initial blowout of 
this bubble into the lower density medium surrounding the 
GMC clump occurs at t ^ 0.03 Myr, and is very aspheri- 
cal due to inhomogeneity of the GMC clump (Fig. [T]). This 
blowout occurs much faster than for a GMC clump with an 
equivalent uniform density - in such a case a ID spherical 
bubble expands to a radius of 4 pc in 0.29 Myr, which is ^ 10 
times slower than for the inhomogeneous GMC clump used 
in our models. Isolated blowouts at a number of distinct po- 
sitions around the surface of the GMC clump rapidly grow 
and merge so that the entire clump is quickly surrounded by 
hot, expanding wind material. As this wind material streams 
out through low density channels, clump material is ablated 
into the flow and the clump gradually loses mass. Fig. [2] 
shows snapshots of the temperature evolution during this 
initial phase at t = 0.03, 0.13 and 0.22 Myr. Here the hot 
wind is clearly visible streaming through the cold GMC. It 
is noticeable that dense parts of the clump can shield and 



protect less dense material in their "shadow", though the 
ability of the hot, high pressure gas to flow around denser 
objects mitigates this effect to some extent. This behaviour 
is apparent in Fig. [l] where material towards the bottom 
right of the clump is protected against ablation from the 
cluster wind by intervening dense material. 

A swept-up shell around the outer edge of the bubble is 
visible in both Figs. ^ and [2] In this simulation the shell is 
reasonably thick and confines the hotter gas in the bubble 
interior, but for bubbles expanding into a higher density 
medium the shell is thinner and less stable. In such situations 
the shell fragments, and mini blowouts occur. These expand 
for a certain distance, before stalling, and merging back into 



the "global" shell (see also [Meaburn, Dyson &: Hartquist 
1988 ). 

At this early stage by far the strongest shock in the 
simulation is the reverse shock near the central stellar clus- 
ter. The reverse shock is not spherical - instead its posi- 
tion is largely defined by the presence and location of the 
dense molecular material nearest to the stellar cluster. The 
shock heated cluster wind experiences a large range of densi- 
ties, temperatures and velocities as it flows out of the GMC 
clump. At times it interacts subsonically with its environ- 
ment, and at other times this interaction is supersonic. Rel- 
atively weaker shocks occur within the shock heated wind 
as it flows into and past the densest molecular gas. The out- 
flow is generally quite turbulent in nature. This can affect 



the rate at which material is stripped from clouds (Pittard 
et al.| '2009). The rate at which clouds lose mass through 
hydrodynamic ablation has been investigated in detail by 
Pittard et al.l (12010). 



3.2 Evolution during the star's MS stage 

Fig. |3] shows the early evolution of the cluster environment 
after the initial blow out has occured but while all three 
stars remain on the MS. The low density channels through 
the GMC clump have left their imprint on the outflowing 
cluster wind as similarly low density channels. These chan- 
nels contain hot (^ 10^ K), fast flowing gas which is rela- 
tively unimpeded by dense gas along its route, with a typical 
velocity of around 1000 kms~^. The flows are mass-loaded 
as denser material along their edges is mixed in. The orien- 
tation of the channels alters slightly as the simulation pro- 
gresses due to the small velocity dispersion few kms~^) 
of the dense molecular material in the GMC clump causing 
its structure to change with time. By about t = 3 Myr the 
position of the channels seems to have settled and they are 
reasonably stable. 

Fig. [4] shows slices of the simulation at t = 0.79 Myr in 
three different planes, while Fig. [5] shows the situation at 
t — 3.41 Myr. The channels carved by the cluster wind are 
also evident in the panels in these figures, and the changes 
in orientation of the channels is apparent between the two 
times. The left panels in Figs. [4] and [5] correspond to the 
xy-plane, which is the default used for the other 2D slices 
within this paper. However, the center and right panels in 
Figs. [4] and |5] show that dense clump material remains closer, 
for longer, to the central star cluster in the xz and yz-planes. 

Together, Figs. [3][5] reveal that the reverse shock ex- 
pands and becomes gradually more spherical with time as 
the cluster wind drives out more material and high density 
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material continues to ablate away. By t — 4Myr the radius 
of the reverse shock has increased to ^ 5 pc, though its ra- 
dius is ~ 3pc at the position of the closest dense cloud to 
the centre of the cluster. 

Fig.[6]shows the pressure at the reverse shock as a func- 
tion of time. The pressure steadily declines throughout the 
main sequence phase of the most massive star, reaching a 
value of 3 X 10~^^dyncm~^ at t = 4Myr. It is difficult to 
know a priori what value to expect for the reverse shock 
pressure. At t = 1 Myr, Eq. 22 of 'Weav er et al.| ( |1977 ) 
gives Pbub = 1.6 X 10"^°dyn cm ^ assuming an ambient 
density of 9 x 10~^^ gcm""^ (roughly the average density of 
our GMC clump), and Pbub = 1-4 x 10~^^dyncm~^ assum- 
ing an ambient density of 3.33 x 10~^^ g cm~^ (the density of 
the medium outside of our GMC clump). Compared to our 
measured pressure of ^ 10~^^ dyncm"^ at this time, we see 
that the former estimate is too high, while the latter is too 



low. Fig. 5 of Harper-Clark & Murray (2009) reveals that 



Pbub is about a factor of 4 — 5 lower, when the covering frac- 
tion Cf ^ 0.3 — 0.6, than the Weaver et al. (1977) estimate. 



This implies that our finite-sized and porous clump has an 
effective covering fraction Cf < 0.3. 

The way that Pbub is set is fundamental to the evolu- 
tion of the flow and its affect on the GMC clump. [Harper-] 
Clark & Murray ( 2009 ) claim that the dynamics of a leaky 



bubble is set by Phii, the pressure of the ionized gas com- 
ponent. Their argument is that when Prs drops to ^ Phii, 
the escape of hot gas through gaps in the bubble shell slows 
as the hot gas is impeded by the cooler gas. However, it 
seems more likely that it is simply determined by the cov- 
ering fraction of the shell Cf and the ram pressure of the 
wind at the shell. If Cf is reasonably high, then individual 
bow shocks around the shell fragments/dense clumps merge 
to create a single reverse shock. The reverse shock will have 
an increasingly large stand-off distance (and thus smaller 
distance from the cluster) as Cf approaches unity. On the 
other hand, if C/ is reasonably low, then the bow shocks 
around individual clumps maintain their identity for longer, 
only merging downstream to create a global reverse shock 
at larger radii from the cluster. Such behaviour can be iden- 



tified in the simulations presented in Pittard et al. (2005) 



and Aluzas et al. (2012) 



3.3 Later evolutionary stages 

3.3.1 Response due to the evolution of the 35 Mq star 

The most massive star evolves to a RSG after 4 Myr. At this 
point its wind speed decreases to Voo = 50kms~^ and its 
mass loss rate increases to M = 10~^ M© yr~^ (see Table [ij. 
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Figure 3. Density slices during the MS phase of the simulation (SimA) in the xy-plane. The last panel shows the density of gas in the 
cluster environment shortly before the most massive star transitions to a RSG. The channels carved by the cluster wind in the GMC 
clump structure slowly evolve over this period. The density scale is shown in the left panel. 
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Figure 6. Pressure at the reverse shock as a function of time. 



This change results in a slower and denser cluster wind. The 
total kinetic power of the cluster wind reduces by about a 
half, from 1.14 x lO^^ergss"^ to 5.87 x 10^^ergss"\ while 
the cluster wind becomes dominated by RSG material. This 
transition is shown in the top rows of Fig.|7|and Fig.js] The 
reverse shock moves inward to reestablish pressure equilib- 
rium with the weaker cluster wind. This depressurises the 
previously shocked gas and leads to a rapid fall in tempera- 
ture of the hottest gas in the simulaton. 

The RSG-enhanced cluster wind is much denser than 
the wind blown when all three stars were on the MS. As it 
interacts with the surrounding gas it is compressed into a 



thin shell which is Ray leigh- Taylor (RT) unstable. RT fin- 
gers are visible in the top right panel in Fig. [7| and Fig. [S] 
These are short lived, lasting approximately 0.04 My r. The 
most massive star remains in the RSG phase for 0.1 Myr, at 
which point the RSG-enhanced cluster wind has expanded 
to a typical radius of ^ 5 pc. 



The most massive star then evolves into a Wolf Rayet 
star, with a mass-loss rate of 2 x 10~^ Mq yr~^ and a wind 
speed of 2000km s~^. This change results in a much faster 
and more powerful cluster wind. The total kinetic power of 
the cluster wind increases by nearly two orders of magnitude 
to 2.59 X 10^^ ergss"^. This transition occurs at t = 4.1 Myr 
and can be seen in the middle row of Fig. [7| and Fig. [S] 
The more powerful cluster wind forcefully pushes back the 
dense RSG material to beyond the position of the reverse 
shock during the previous MS phase. The typical radius of 
the reverse shock increases from about 5pc at t = 4.14 Myr 
to ~ 8pc at t = 4.4 Myr (see middle row of Fig. [7|. The 
shocked cluster wind is ~ 10^ times hotter than was the 
case when the cluster wind was "RSG-enhanced". Hot gas 
pervades almost completely the computational volume by 
t = 4.15 Myr (see Fig.[8|. 
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3.3.2 Impact of the first supernova explosion 

After 4.4 Myrs of evolution the most massive star explodes 
as a supernova (see the bottom rows of Figs. [t] and |8|, adding 
10 M0 of ejecta and 10^^ ergs of energy into the surround- 
ings. The SN ejecta sweeps up the WR-dominated cluster 
wind into a thick shell, which propagates at high speed 
through the lower density regions surrounding the explo- 
sion, heating the gas to very high temperatures. The for- 
ward shock has already propagated off the grid by the time 
of the snapshot in the bottom left panel of Fig.|7| The highly 
aspherical reverse shock seen in this panel is caused by the 
forward shock reflecting off dense clumps of gas in the sim- 
ulation. The reverse shock then moves inwards towards the 
central cluster and nearly overwhelms the MS-dominated 
cluster wind from the remaining two stars (bottom middle 
panel of Fig.[7|, which is seen sweeping up the cold SN ejecta 
in the bottom left panel. The reverse shock of the cluster 
wind is forced back to a radius of ~ 0.9 pc by t = 4.42 Myr, 
but slowly begins to expand to ^ 1.5 pc by t = 4.47 Myr 
as the pressure within the supernova remnant diminishes. 
One of the main effects of the SN is to fill the low-density 
wind-carved channels through the GMC clump with denser 
material which the cluster wind must again clear. The very 
densest fragments left over from the original GMC clump 
are not only resistant to ablation from the cluster wind, but 
also are affected very little by the propagation of the SN 
Shockwave. 

Fig. |9^) shows the temperature of the gas in the sim- 
ulation during the 5000 years after the most massive star 
goes supernova. The green dashed line shows the state just 
before the supernova explosion, and the thick, red, solid line 
corresponds to just afterwards. The maximum temperature 
on the grid jumps by almost an order of magnitude in the 
first 100 years (from ^ 10^ to ^ 4 x 10^ K) and the amount 
of material between 10^ — 10^ K increases from about 1 M© 
to 10 M0. The maximum temperature continues to increase 
up to 10^ K as the shockwave propagates through the cloud. 
Gas with T < 100 K responds more slowly to the SN explo- 
sion, but it is clear that some of it is heated to T ^ 10^ K 
by the passage of the shock wave. A significant proportion 
of the coldest gas with T < lOK is hardly affected, how- 
ever. This very cold gas is situated in the very densest re- 
gions of the remaining GMC clump. As the shockwave passes 



through the simulation, the less dense material surrounding 
these regions is heated and ionized, but the densest parts 
are relatively untouched. This is because the transmitted 
shock speed through the densest clumps is as low as a few 
lO's of kms~^, and the gas cools back to its original pre- 
shock temperature on timescales as short as a few 10 's of 
yrs. More typical cooling times are 10^ yrs or so, but it is 
clear that the transmitted shocks into the densest clumps 
are highly radiative (the crossing time of the transmitted 
shocks through these fragments is ^ 10^ yrs). The strong 
cooling of this gas is responsible for the rapid rise in the 
H2 mass following its initial drop after each SN explosion 



(see Sec. 3.6). The SN ejecta first begins to leave the grid 
approximately 4000 years after the explosion. 

Fig.|9]3) shows the pressure (p/k) during this time. Af- 
ter the SN explosion the pressure increases by four orders of 
magnitude from 10^ - 10^^ K cm . The maximum pressure 
then slowly decreases as the remnant expands adiabatically. 
By t — 4.4046 Myr (black double-dashed line) the maximum 
pressure on the grid is back to ~ 10^ K cm"*^, but now there 
is ^ 860 M© of material at lO'^ ^ p/k ^ lO^Kcm"^, com- 
pared to the 0.02 M© of gas in this pressure range prior to 
the SN. 50, 000 yrs after the SN explosion the net effect is a 
shift in gas from lower (10^ < p/k < 10^ Kcm~^) to higher 
(10^ < p/k < 10^ Kcm"^) pressure. 

The maximum gas velocity also increases rapidly in 
the first 100 years after the supernova explosion, as shown 
m Fig. jg]:). The maximum gas velocity increases from 
2000 km s~^ to 10, 000 km s~^ as the hot, high pressure ejecta 
starts its expansion. Over the next 5000 yrs the maximum 
velocity drops to ^ 5000 km s~^ and there is significantly 
more mass with v ^ lOOkms"^ than was the case pre- 
SN. The majority of the gas continues to have a velocity of 
1 — 10 km s~^ , however, which highlights the weak coupling 
of the SN explosion to the densest parts of the surrounding 
gas. 



3.3.3 Response due to the evolution of the 32 Mq and 
28 Mq stars 

0.1 Myr after the explosion of the most massive star, the 
32 Mq star evolves off the MS and onto the RSG branch (see 
Table [T]) , decreasing the kinetic power of the cluster wind 
still further to 2.7 x 10^^ ergss"^. This star follows the same 
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Figure 7. Density slices from SimA in the xy-plane during the RSG, WR and SN stages of the highest mass star. This star transitions 
from the MS to the RSG stage at t = 4.0 Myr (top middle panel), from the RSG to the WR stage at t = 4.1 Myr (top right panel), and 
explodes at t = 4.4 Myr (bottom left panel). The other two stars remain on the MS during this time. The density scale is shown in the 
top left panel. 



post-MS evolution as the most massive star, and blows a 
dense slow moving RSG wind followed by a high momentum 
WR wind before finally exploding. The SN explosion again 
imparts 10 of material and 10^^ ergs s~^ of energy into 
the surroundings. After a further 0.1 Myr the lowest mass 
(and only remaining) star considered in our model evolves 
onto the RSG branch, reducing the kinetic power of the 
cluster wind to 7.9 x 10^^ ergs s~^. This star explodes as 
a supernova at t = 5.4 Myr. At this moment there are no 
wind-blowing stars remaining in our model, and the ambient 
and enriched gas gradually disperses and flows off the grid 
boundaries. 



3.4 Mass and energy fluxes into the wider 
environment 

Fig. [To] shows the total mass flux off the grid as a func- 
tion of time. The mass flux is zero until the first blowouts 
reach the edge of the grid and then steadily increases up 
until t ^ 0.4 Myr as various shells of swept up mate- 
rial reach the outer boundary of the simulation. The ma- 
terial behind the shells is less dense, and therefore the 
mass flux declines as the shells leave the grid. The mass 
flux stabilizes at t ~ 1.15 Myr and then begins a slow, 
almost linear, increase from 1.6 x lO~^M0yr~^, reaching 
nearly 3 x lO~^M0yr~^ at t = 4.0 Myr. In comparison. 
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Figure 8. Temperature slices from SimA in the xy-plane during the RSG, WR and SN stages for the highest mass star. The first panel 
corresponds to just after this star has become a RSG. The star transitions to a WR star ai t = 4.1 Myr (left panel middle row) and 
explodes as a SN at t = 4.4 Myr (left panel bottom row). The other stars remain on the MS during this time. The temperature scale is 
shown in the left panels. 



the mass-loss rates of the three stars during the MS is 
5 X 10"'^, 2.5 X 10"'^ and 1.5 x IQ-'^Moyr"^ for the 35 M©, 
32 M0 and 28 Mq stars, respectively, giving a cluster mass- 
loss rate of 9 x 10~^ yr~^. This indicates that the cluster 
wind is "mass-loaded"[jby factors of 200 — 300 as it streams 
through the molecular material in the GMC clump, and that 



^ Although we use the term "mass-loaded" here, the molecular 
material ablated by the cluster wind is not fully mixed into the 
fiow by the time that it leaves the grid. Since distinct phases are 
still identifiable, "mass entrainment" may be a more appropriate 
term. 



less than 1 per cent of the material leaving the grid during 
this period originated in the stellar winds. 

There is a slight decrease in the mass flux off the grid be- 
tween t = 4.0 and 4.1 Myr, during the RSG stage of the high- 
est mass star in which the flow de-pressurizes. This is fol- 
lowed by an increase in the mass flux to ~ 7 x 10~^ M© yr~^ 
during the subsequent WR stage. This represents a "mass- 
loading" factor 35 X the mass- loss rate of the cluster wind 
during this period (2.04 x 10~^ M© yr~^). Therefore, only 3 
per cent of the material leaving the grid is from the stellar 
winds. The mass flux jumps to 2.5 x 10"*^ M© yr~^ following 
the first SN explosion, and peaks at 3.9 x 10""^ Mq yr~^ fol- 
lowing each of the second and third explosions. In between 
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Figure 9. Histogram of the gas temperature (top), pressure (mid- 
dle), and velocity (bottom) during the first supernova explosion. 
In each panel the solid red line corresponds to a time just after 
the explosion occurs. 



the SN explosions, the mass flux peaks at 2.0 x 10~^ yr~^ 
at t = 5.31 Myr. Clearly the WR winds and SN explosions 
act to speed up the rate at which molecular material is 
cleared from the cluster surroundings. 

The top panel of Fig. [iT] shows the total energy flux 
off the grid as a function of time. The energy flux increases 
from zero at the time when the shell first encounters the 
grid boundary, to ^ 8.5 x 10^^ ergss"^ at t = 1 Myr, with a 
gradual and linear decline to 7.5 x 10^^ ergs s~^ at t = 4 Myr. 
The linear decline correlates with the gradual increase in the 
ablation rate during this period. These values for the energy 
flux compare to the total kinetic power of the cluster wind 
of 1.14 X 10"^^ ergss"^. It is clear, therefore, that the shocked 
cluster wind is largely adiabatic, but that nevertheless about 
one quarter to one third of the injected wind power is lost to 
radiative processes. It is interesting to note that |Bruhweiler| 
et al. (2010) invoke substantial radiative losses due to the 



turbulent interaction of stellar winds with inhomogeneous 
surroundings in order to explain their observations of the 
Rosette Nebula. 

The combined energy input by the 3 0-stars during the 
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Figure 10. Total mass flux off the grid as a function of time. 

first 4 Myr is 1.44 x 10^° ergs. When the most massive star 
evolves to the RSG phase there is a decrease in the kinetic 
power of the cluster wind to 5.87 x 10^^ ergs s~^, and this 
is reflected in a corresponding decrease in the energy flux 
off the grid, which, however, drops below this value. This 
"overshoot" is likely due to some combination of enhanced 
cooling in the denser cluster wind and possible overstable 
behaviour of the flow. The energy flux off the grid increases 
by two orders of magnitude to 1.8 x 10^^ ergs s~^ when the 
most massive star has evolved to a WR star. Comparing to 
the kinetic power of the cluster wind at this time (2.59 x 
10^^ ergs s~^), we find that radiative losses are again about 
30 per cent. The combined energy input by the WR star and 
the 2 remaining 0-stars during the period 4.1 — 4.4 Myr is 
2.46 X 10^° ergs. Altogether, the most massive star injects 
^ 3.9 X 10^° ergs of energy via its wind during its lifetime, 
and about 70 per cent of this escapes to large distances. 

The bottom panel of Fig. [TT] shows the total energy 
flux off the grid around the time of the first SN explosion 
which occurs at t = 4.4 Myr. The energy flux rises steeply 
as the blast shock propagates off the grid, and peaks at ^ 
5 X 10^^ ergs s~^ about 5000 yr after the explosion. Thereafter 
the energy flux steadily decreases. The integrated energy off 
the grid during this time reveals that greater than 99 per 
cent of the SN energy flows off the grid, and less than 0.5 per 
cent is radiated. Clearly the SN energy propagates through 
the environment relatively unimpeded by the dense clumps. 

The subsequent energy flux from the grid is dominated 
by the WR stages from the two lower mass stars and their 
subsequent explosions, and drops below 10^^ ergss"^ as relic 
hot gas expands and dissipates. 



3.5 Evolution of column densities 

Fig. [12] shows the time evolution of the average column den- 
sity, iVn, from the centre of the cluster. This is calculated 
over 10^ individual sight lines spaced equally in solid angle 
and traced out to the edge of the grid. The column den- 
sity is greatest at t = 0, but diminishes with time as the 
stellar winds push the clump material away from the clus- 
ter (some of this decrease is also due to material leaving 
the grid). A factor of 100 reduction from an initial value of 
4 X 10^^ cm~^ occurs by 4 Myr. Then, as the most mas- 
sive star enters its RSG stage, iVn increases by over 1 dex 
to nearly lO^^cm"^. This increase is short-lived however, 
because the density of the cluster wind drops significantly 
once the most massive star enters its WR stage. Immediately 
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Figure 11. Total energy flux off the grid as a function of time. 
Top: Entire simulation. Bottom: Focussed on the flrst SN ex- 
plosion. The steps are caused by the cadence of the timesteps 
analysed. 




Figure 12. The evolution of the average column density from 
the centre of the cluster. 



prior to the first SN explosion, Nu ^ 2 x lO^^cm"^. The 
10 of ejecta from the explosion momentarily increases 
the average column density to more than lO^^cm"^, but 
this rise is extremely short-lived and less than a few hun- 
dred years in duration. Nu then gradually increases as the 
cluster wind again begins to fill the nearby environment with 
mass, especially once the second most massive star enters its 
RSG stage. Similar behaviour in the evolution of Nh then 
occurs as the remaining massive stars evolve in turn through 
their various wind and SN stages. 

Fig. [13] shows histograms of the column density dis- 
tribution at specific times in the simulation. At t = Oyr, 
columns up to ~ 10^^ cm~^ exist for sightlines through the 
densest parts of the initial GMC clump, but a very small 
fraction of sightlines from the cluster experience column 
densities as low as 10^°cm~^. At t = 3.99 Myr, the dis- 
tribution of column densities is much broader, with those 
passing through dense and relatively nearby regions having 



teresting to see how the distribution has changed shape by 
t — 4.02 Myr, once the most massive star has entered its 
RSG stage. The RSG-enhanced cluster wind "fills in" all of 
the low column density sightlines so that the minimum value 
is now A^H ~ 3 X 10^° cm~^. The impact of the WR-enhanced 
cluster wind is seen oX t — 4.39 Myr. The higher speed and 
reduced density of the cluster wind now reduces the low 
column density sightlines to ~ 10^^ cm~^. The column den- 
sity distribution at a time just after the first SN explosion is 
shown at t = 4.40000 Myr. The relatively dense ejecta causes 
the minimum value of A^h to rise to ~ 10^^ cm~^, yet this 
drops to just 10^^ cm~^ 1100 yrs later. 

Fig. [M] shows Hammer projections of the column den- 
sity at specific times. The directions with the highest (low- 
est) initial column density maintain their positions through- 
out the simulation. It is clear that the regions of highest col- 
umn density become increasingly "porous" and "shredded" 
by the action of the cluster wind and supernovae. The "fill- 
ing in" of the column density during the time of the RSG- 
enhanced cluster wind (t = 4.02 Myr in Fig. 14) is readily 
apparent. 

Though we do not include photoionization in this work, 
it is interesting to estimate the fraction of ionizing photons 
which could escape to large distances. Along each sight line 
we evaluate whether the integral of A^nr^ aBTiidr exceeds Sd^ 
where r is the radial distance to the stellar cluster, rie is the 
electron number density, is the case B recombination co- 
efficient and Sc\ is the rate of ionizing photons from the stel- 
lar cluster. If this condition is satisfied the ionization front is 
trapped on the grid in that direction. Fig. |15| shows the re- 
sult of this calculation. Initially the ionization front is almost 
completely trapped within the GMC clump surrounding the 
stellar cluster, with less than 1 per cent of ionizing pho- 
tons escaping (as indicated by the red colour). However, by 
t = 0.95 Myr, the stellar winds have sufficiently cleared away 
dense molecular and atomic material that about 40 per cent 
of the ionizing radiation escapes. This increases to nearly 60 
per cent by t = 3.99 Myr. The RSG-enhanced cluster wind 
is then sufficiently dense to completely prevent any ionizing 
radiation escaping to large distances. Once the lower density 
WR-enhanced cluster wind clears the RSG dominated ma- 
terial off the grid the escape fraction increases once more, 
reaching nearly 65 per cent just prior to the first SN explo- 
sion and 75 per cent at t = 4.94 Myr. 



Recently, Pellegrini et al. (2012) determined that the 



column densities up to 10^^ cm ^ and those passing through 
the lowest density material having A^h ~ 10^^ cm~^. It is in- 



mean, luminosity weighted, escape fraction of ionizing pho- 
tons from the HII region population in the LMC and SMC 
is > 0.42 and > 0.40, respectively. Our results are consis- 
tent with these values, and similar values from |Reines et ah] 
( [2008| ,but a more rigorous and detailed analysis of the pho- 
toionization beyond this work is necessary. 



3.6 Comparison of SimA and SimB and evolution 
of the molecular mass 

Fig. [16] shows the initial destruction of the GMC clump for 
SimA and SimB. The higher density clump and increased 
clump radius in SimB inevitably delays the blowout and sub- 
sequent expansion of the cluster wind. However, the nature 
of the interaction is broadly similar, and the same general 
features are seen as the cluster wind percolates through the 
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Figure 13. The evolution of the average column density from the centre of the cluster. From left to right and top to bottom the plots 
are at times t = 0, t = 3.99, t = 4.02, t = 4.39, t = 4.40000 and t = 4.40110 Myr. The numerical value in each panel notes the value of 
logio Nh at which 50 per cent of the sky has a smaller or larger column. 




Figure 14. Hammer projections of the column density at specific times. 
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"porous" molecular environment, including the formation of 
fast-flowing low-density channels. 

The mass flux off the grid from SimB behaves qual- 
itatively similarly to that from Sim A, although the exact 
values are a little higher. There is again a broad peak 
caused by the initial encounter of the shell with the bound- 
aries, in which the mass flux off the grid peaks at a rate 
of 4.6 X lO"^M0yr"^ at t = 0.67 Myr. The mass flux de- 
clines to a minimum of 2.9 x 10~^Mq yr~^ at t ~ 1.6 Myr, 
and then increases roughly linearly with time to reach a 
value of 7 X 10"^ M© yr~^ at t = 4 Myr. The latter indicates 
a "mass-loading" factor of nearly 1000. The mass flux in- 
creases to 1.6 X 10~^ M0 yr~^ during the first WR stage, 
a mass- loading factor of nearly 80. A peak mass flux of 

4.3 X lO"^M0yr~^ is attained following the first SN ex- 
plosion, with peaks at 7 and 8 x 10~^ M© yr~^ following the 
second and third explosions respectively. 

The energy flux off the grid from SimB peaks at ~ 

7.4 X 10^^ ergs s~^ at t = 1.5 Myr, followed by a gradual 
decline to a roughly constant rate of ^ 6.4 x 10^^ ergs s~^ 
between t = 2.5 — 4 Myr. During this latter period, nearly 
half of the cluster wind power is radiated. The energy flux 
plateau's at about 1.45 x 10^^ ergs s~^ during the first WR 
stage (indicating that again nearly half of the input energy is 
being radiated) , but similar values are reached only towards 
the very end of the other two WR stages. 

The evolution of the total H2 mass throughout the simu- 
lation volume for both SimA and SimB is shown in Fig. |17| 
Due to the larger clump radius and higher density, SimB 
has approximately three times the initial H2 mass as SimA. 
At very early times during the initial blowout the H2 mass 
decreases by ^ 10% as molecular gas at relatively low densi- 
ties is cleared out of the GMC clump. Both simulations then 
show the same general trend of steady H2 depletion during 
the MS phase of the cluster wind. SimA loses H2 mass at a 
rate of 1.74 x lO"^M0yr~^ between t = 1 - 4Myr. Since 
this is comparable to the mass flux off the grid during this 
time it further reinforces the point that most of the mass 
streaming into the wider environment was orginally molec- 
ular material. 

The evolution of the H2 mass becomes more exciting 
when the stars undergo evolutionary transitions. There are 
slight rises during periods when the cluster wind is dom- 
inated by RSG mass-loss, and more significant decreases 
during WR dominated periods when the molecular mate- 
rial is ablated at a much faster rate. But most dramatic of 
all is the response of the H2 gas to a supernova explosion. 
Immediately after an explosion, the H2 mass undergoes a 
rapid and steep decline which is quickly followed by a recov- 
ery to a similar or higher value than the pre-SN H2 mass. 
The depletion of H2 at this time is caused by its conversion 
to neutral and ionized hydrogen as it is heated by the SN 
Shockwave which propagates relatively slowly through such 
dense molecular material. Note that it is not due to the H2 
gas being expelled from the grid. Once the shockwave has 
passed the densest shocked (neutral and atomic) gas cools 
and reverts quickly to its previous molecular state, result- 
ing in the replenishment of the H2 mass seen in Fig. |17| 
The overshoot of the H2 mass relative to its pre-SN value is 
likely a result of the higher gas pressure which now exists 
within the cluster environment, and which allows some of 
the originally atomic hydrodgen to become molecular. The 



mass of molecular gas resumes its decline once the next most 
massive star in the cluster enters its WR stage. 

Since most of each SN's energy escapes along the low 
density channels through the remains of the GMC clump 



(see Sec. 3.4), the SNe couple very weakly to the densest. 



low volume filling factor, gas. Hence the very densest regions 
are relatively immune to the effects of the supernovae (see 
also Figs, [t] and [8|. Indeed, Fig. 17 shows that overall, the 
supernova shocks actually increase the amount of molecular 
material in the cluster environment. This rather unexpected 
result demonstrates that the inhomogeneity of the environ- 
ment, the density range of the gas, and its effective porosity, 
are key considerations in understanding massive star feed- 
back. This picture also differs substantially from the many 
models of mass-loaded supernova remnants in the literature 
which envisage the complete mixing of material injected into 



the remant from embedded clumps (e.g. White & Long 1991 
Dyson et al.||2002| [Pittard et al.|2003t . 

The last supernova occurs at t = 5.4 Myr, after which 
there is ~ 470 Mq (14.4% of the original mass) of H2 remain- 
ing in SimA and ~ 3100 M© (29.8% of the original mass) in 
SimB. 

Fig.|18|shows the mass of H2 contained within the initial 
clump radius (4 and 5pc for SimA and SimB respectively). 
The depletion of H2 is much more rapid within this radius as 
the cluster wind not only ablates but also pushes molecular 
material away from the stars. By the time the lowest mass 
star explodes there is very little H2 present within the orig- 
inal volume of the GMC clump in either simulation (< 0.1 
per cent of the original mass for SimA). The features due 
to the evolution of the stars shown in Fig. [T7| are apparent 
here also, although to a lesser degree. 



4 CONCLUSION 

This paper investigates the effects of massive star feedback, 
via stellar winds and supernovae, on the inhomogeneous 
molecular environment left over from the formation of a stel- 
lar cluster from a GMC clump. The remains of the GMC 
clump confines and shapes the initial structure of the ex- 
panding wind-blown bubble, which breaks out of the clump 
along paths of least resistance. Hot, high speed gas flows 
away from the cluster through low-density channels opened 
up by these blow-outs. Mass is loaded into these flows from 
the ablation of dense clumps embedded within them and 
from material stripped from the dense gas which confines 
and directs the flows. This complex interaction of the clus- 
ter wind with its environment is far removed from the results 
of simple spherically symmetric models. Increasing the den- 
sity (by a factor of two) and the radius (by 25 per cent) of 
the cluster-forming GMC clump does not significantly affect 
this qualitative picture. 

The density, temperature, pressure and velocity of gas 
in the cluster environment all span many orders of magni- 
tude. The hottest gas typically occurs at the reverse shock 
of the cluster wind, and cools as it expands away from the 
cluster and mixes in with denser surrounding material. A 
multitude of weaker shocks exist around dense inhomogeni- 
ties entrained into gas flowing at mildly supersonic speeds. 
In addition, gas within the cluster environment is subject 
to changes of several orders of magnitude in the dynamic 
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Figure 16. Comparison of the initial expansion of the cluster 
wind in the xy-plane. [Left]: Sim A and [Right]: Sim B. Sim B 
has a higher ambient density and a larger clump radius than Sim 
A. The density scales for each simulation are shown in the top 
panels. 



pressure as the stars in the cluster evolve through their MS, 
RSG and WR stages, and explode as supernovae. 

Our simulations show how molecular material is gradu- 
ally ablated by the cluster wind and pushed away from the 
stellar cluster. Despite mass-loading or entrainment factors 
of several hundreds during the MS phase of the cluster wind, 
and several tens during the later WR-dominated phase, the 
destruction and sweeping up of molecular material is a rel- 
atively slow process, and a substantial amount of molecular 
mass remains when the first star explodes as a SN. We find 
that the shocks resulting from SN explosions couple very 
weakly to the molecular gas, due to its small volume fill- 
ing fraction, and the ease with which the energy from the 
SN can "by-pass" it. The high porosity of the GMC clump 
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Figure 18. Mass of H2 contained within the initial GMC clump 
radius of SimA (4pc, solid red line) and SimB (5pc, dotted green 
line). 

at this stage allows the SN blast to rip through the clus- 
ter in a largely unimpeded fashion, with the forward shock 
refracting around dense inhomogeneities. The early evolu- 
tion of the remnant is markedly different from the standard 
spherically symmetric picture. Although the SN shock de- 
stroys molecular material which it overruns, we find that the 
cooling times of the densest regions are very short and allow 
molecular material to quickly reform. At least in our simula- 
tions, the stellar winds appear to be a more effective agent at 
removing molecular material from the cluster environment, 
despite injecting less energy than the SNe. 

Examination of the energy flux off the hydrodynamic 
grid reveals that between one quarter and one half of the 
energy injected by the stellar winds is radiated away, with 
the remaining energy available to do work on the immediate 
surroundings. In comparison, more than 99 per cent of the 
energy from each SN explosion escapes into the wider envi- 
ronment. These fractions are clearly dependent on the initial 
conditions of our models. In particular, our limited study 
suggests that the fractions of energy radiated away will in- 
crease for denser cluster environments, and vice- versa. We 
are performing additional calculations to determine the con- 
ditions necessary for almost all of the injected energy to be 
radiated. This is relevant to some of the super star clusters 
in M82 where the thermalization effiency of the stellar feed- 
back (winds plus supernovae) does not exceed a few percent 
( jSilich et al.'|2QQ9t . 

We also estimate the fraction of ionizing photons which 
escape the cluster environment. This increases from less than 
1 per cent at the start of the simulation, to 40 per cent when 
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the cluster is 1 Myr old, and to 60 per cent after 4 Myr. The 
escape fraction momentarily reduces when the cluster wind 
is RSG dominated, but the overall trend is of an increasing 
escape fraction as dense material is pushed away from the 
cluster. 

In a future paper we will examine synthetic emission cal- 
culated from our simulations which we will compare against 
observations. We will also address some of the simplifica- 
tions of our current model such as the neglect of radiative 
pressure and photoionization. 
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